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ABSTRACT 

We present results of large-scale three-dimensional simulations of supersonic Euler turbulence with the piece- 
wise parabolic method and multiple grid resolutions up to 2048 3 points. Our numerical experiments describe 
non-magnetized driven turbulent flows with an isothermal equation of state and an rms Mach number of 6. We 
discuss numerical resolution issues and demonstrate convergence, in a statistical sense, of the inertial range 
dynamics in simulations on grids larger than 512 3 points. The simulations allowed us to measure the absolute 
velocity scaling exponents for the first time. The inertial range velocity scaling in this strongly compress- 
ible regime deviates substantially from the incompressible Kolmogorov laws. The slope of the velocity power 
spectrum, for instance, is -1.95 compared to -5/3 in the incompressible case. The exponent of the third-order 
velocity structure function is 1 .28, while in incompressible turbulence it is known to be unity. We propose a nat- 
ural extension of Kolmogorov's phenomenology that takes into account compressibility by mixing the velocity 
and density statistics and preserves the Kolmogorov scaling of the power spectrum and structure functions of 
the density-weighted velocity v = p 1 ' 3 !!. The low-order statistics of v appear to be invariant with respect to 
changes in the Mach number. For instance, at Mach 6 the slope of the power spectrum of v is -1.69, and the 
exponent of the third-order structure function of v is unity. We also directly measure the mass dimension of 
the "fractal" density distribution in the inertial subrange, D m « 2.4, which is similar to the observed fractal 
dimension of molecular clouds and agrees well with the cascade phenomenology. 

Subject headings: hydrodynamics — instabilities — ISM: structure — methods: numerical — turbulence 



1. INTRODUCTION 

Understanding the nature of supersonic turbulence is of 
fundamental importance in both astrophysics and aeronau- 
tical engineering. In the interstellar medium (ISM), highly 
compressible turbulence is be lieved to control star form ation 
in dense molecular clouds dPadoan & Nord lund 2002J). In 
radiation-driven outflows from carbon-dominant Wolf-Rayet 
stars, supersonic turbulence creates highly clumpy structure 
that is stochastica lly variable on a very short time-scale (e.g., 
lAcker et al. 2002). Finally, a whole class of more terrestrial 
applications deals with the drag and stability of projectiles 
traveling through the air at hypersonic speeds. 

Molecular clouds have an extremely inhomogeneous struc- 
ture and the intensity of their inter nal mot i ons co rresponds to 
an rms Mach number of order 20. iLarsonl (119811) has demon- 
strated that within the range of scales from 0.1 to 100 pc, the 
gas density and the velocity dispersion tightly correlate with 
the cloud size. 2 Supported by other independent observational 
facts indicating scale invariance, these relationships are often 
interpreted in terms of supersonic turbu lence with character - 
istic Reynolds numbers Re~ 10 8 (Elmegree n & Scaloll2004l) . 
Within a wide range of densities above 10 3 cm -3 , the gas tem- 
perature remains close to as 10 K, since the thermal equili- 
bration time at these densities is shorter than a typical hy- 
drodynamic (HD) timescale. Thus, an isothermal equation of 
state can be used as a reasonable approximation. Self-gravity, 
magnetic fields, chemistry, cooling, and heating, as well as 
radiative transfer, should ultimately be accounted for in turbu- 
lent models of molecular clouds. However, since highly com- 

1 Also at: Sobolev Astronomical Institute, St. Petersburg State University, 
St. Petersburg, Russia. 

2 For an earlie r versi on of what is now k nown as Larson's relations, see 
IKaplan & Pronikl 119531) and also see lBrunfl 120031) for the latest results on 
the velocity dispersion-cloud size relation. 



pressible turbulence still remains an unsolved problem even in 
the absence of magnetic fields, our main focus here is specif- 
ically on the more tractable HD aspects of the problem. 

Magnetic fields are important for the general ISM dynam- 
ics and, particularly, for the star formation process. Ob- 
servations of molecular clouds are consistent with the pres- 
ence of super- Alfvenic turbulence (Padoan et al. 2004a|), and 
thus, the average magnetic field strength may be much 
smaller than require d to support the clouds against t he grav- 
itational collapse dPadoan & Nordlund! [19971 1 19991) . Even 
weak fields, however, have the potential to modify the proper- 
ties of supersonic turbulent flows through the effects of mag- 
netic tension and magnetic pressure and introduce small-scale 
anisotropics. Magnetic tension tends to stabilize hydrody- 
namically unstable post-shock shear layers dMiura & Pritchettl 
119821: iKeppens et all 19991: iRvu et al.ll2000t lBatv et alj|2003l) . 
The shock jump conditions modified by magnetic pressure re- 
sult in substantially different predictions for the initial mass 
function of stars forming in non- magnetic and magn etized 
cases via turbulent fragmentation (Pad oan et al.ll2007|). Re- 
mark ably, although not surprisingly (e.g.. lArmi & Flamenj 
1985), the impact of magnetic field on the low-order statistics 
of super- Alfvenic turbulence appears to be rather limited. At a 
grid resolution of 1024 3 points, the slopes of the power spec- 
tra in our simulations show stronger sensitivity to the numer- 
ical diffusivity of the scheme of choice t han to the presence 
of the magnetic field dPadoan et al.ll2007h . The similarity of 
non-magnetized and weakly magnetized turbulence will allow 
us to compare our HD results with those from super- Alfvenic 
magnetohydrodynamic (MHD) simulations where equivalent 
pure HD simulations are not available in the literature. 

Numerical simulations of decaying supersonic HD tur- 
bulence with the piecewise parabolic method (PPM) in 



2 



KRITSUK, NORMAN, PADOAN, & WAGNER 



5 




7 8 9 10 

t/t d 

FIG. 1 . — Time evolution of the rms Mach number in the 1024 3 simulation 
of driven Mach 6 turbulence. 



two dimensions were pioneered by Pas sot et all {l988) 3 
and then followed up with hi gh-resolutio n two- a nd three- 
dimensional simulations by iPorter et ail dl992alfbl 11994 
119981) . ISvtine et all (120001) compared the results of PPM Eu- 
ler computations with PPM Navier-Stokes results and showed 
that the Euler simulations agree well with the high-Re limit at- 
tained in the Navier-Stokes models. The convergence in a sta- 
tistical sense as well as the direct comparison of structures in 
configuration space indicate the ability of PPM to accurately 
simulate turbulent flows ov er a wide range of scales. More 
recently, IPorter et"aT] d2002l) discussed measures of intermit- 
tency in simulated driven transonic flows at Ma ch numbers 
of the order unity on grids of up to 512 3 points. Por ter et all 
(1999) review the results of these numerical studies, focusing 
on the origin and evolution of turbulent structures in physical 
space as well as on scaling laws for two-point structure func- 
tions. One of the important results of this fundamental work is 
the demonstration of the compatibility of a Kolmogorov-type 
(iKolmogorovl 1 1 94 TaL K41) spectrum with a mild gas com- 
pressibility at transonic Mach numbers. 

Since most of the computations discussed above assume a 
perfect gas equation of state with the ratios of specific heats 
7 = 7/5 or 5/3 and Mach numbers generally below 2, the 
question remains whether this result will still hold for near 
isothermal conditions and hypersonic Mach numbers charac- 
teristic of dense parts of star-forming molecular clouds where 
the gas compressibility is much higher. What kind of coher- 
ent structures should one expect to see in highly supersonic 
turbulence? Do low-order statistics of turbulence follow the 
K41 predictions closely in this regime? How can the statis- 
tical diagnostics traditionally used in studies of incompress- 
ible turbulence be extended to reconstruct the energy cascade 
properties in supersonic flows? How can we measure the in- 
trinsic intermittency of supersonic turbulence? Many of these 
and similar questions can only be addressed with numerical 
simulations of sufficiently high resolution. The interpretation 
of astronomical data from new surveys of the cold ISM and 
dust in the Milky Way by the Spitzer and Herschel Space Ob- 
servatory satellites requires more detailed knowledge of these 
basic properties of supersonic turbulence. 

In this paper we report the results from large-scale numeri- 
cal simulations of driven supersonic isothermal turbulence at 
Mach 6 with PPM and grid resolutions up to 2048 3 points. 
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FIG. 2. — Same as Fig.rjJ but for the maximum gas density. 



The paper is organized as follows: Section [2] contains the de- 
tails of the simulations' setup and describes the input param- 
eters. The statistical diagnostics, including power spectra of 
the velocity, the kinetic energy and the density, and veloc- 
ity structure functions, together with a discussion of turbulent 
structures and their fractal dimension are presented in Section 
[3] In Section [3791 we combine the scaling laws determined 
in our numerical experimen ts to verify a simple compressible 
cascade model proposed by iFleckl (119961) . We also introduce 
a new variable v = p l ^u that controls the energy transfer rate 
through the compressible cascade. We then summarize the 
results and discuss possible ways to validate our numerical 
model with astrophysical observations in Sections @]and|5] 

2. METHODOLOGY 

We use PPM implemented in the Enzo code 4 to solve the 
Euler equations for the gas density p and the velocity u with a 
constant external acceleration term F(x) such that (F(x)) = 0, 



d,p+v-( P u) = o, 

d,u+u-Vu = -Vp/p+F, 



(1) 
(2) 



3 See a review on compressible turbulence by Pouquet et al] 1199 ID for 
references to earlier works. 



in a periodic box of linear size L = 1 starting with an initially 
uniform density distribution p(jc, t = 0) = po(x) = 1 and assum- 
ing the sound speed c = 1 . The equations that were actually 
numerically integrated were written in a conventional form 
of conservation laws for the mass, momentum, and total en- 
ergy that is less compact, but nearly equivalent to equations 
(fl} and (0, since in practice we mimic the isothermal equa- 
tion of state by setting the specific heats ratio in the ideal gas 
equation of state very close to unity, 7=1 .001 . 

The simulations were initialized on grids of 256 3 or 512 3 
points with a random velocity field u(x, t = 0) = uq(x) oc F that 
contains only large-scale power within the range of wavenum- 
bers k/k m j„ £ [1,2], where k m i„ = 2ir, and that corresponds to 
the rms Mach number M. = 6. The dynamical time is hereafter 
defined as t d = L/(2M). 

2.1. Uniform Grid Simulation at 1024 3 

Our major production run is performed on a grid of 1024 3 
points that allowed us to resolve a portion of the uncontami- 
nated inertial range sufficient to get a first approximation for 
the low-order scaling exponents. 

4 See http://lca.ucsd.edu/ 
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We started the simulation at a lower resolution of 512 3 
points and evolved the flow from the initial conditions for five 
dynamical times to stir up the gas in the box. Then we dou- 
bled the resolution and evolved the simulation for another 5td 
on a grid of 1024 3 points. 

The time-average statistics were computed using 170 snap- 
shots evenly spaced in time over the final segment of 4 td- 
(We allowed one dynamical time for flow relaxation at high 
resolution, so that it could reach a statistical steady state after 
regridding.) We used the full set of 170 snapshots to derive the 
density statistics, since the density field displays a very high 
degree of intermittency. This gave us a very large statistical 
sample, e.g. 2 x 10 11 measurements were available to deter- 
mine the probability density function (PDF) of the gas density 
discussed in Section [3~2l The time-average power spectra dis- 
cussed in Sections 13.31 13.51 13.61 and 13.91 are also based on 
the full data set. The velocity structure functions presented in 
Section [3~4l are derived from a sample of 20% of the snapshots 
covering the same period of 4 td- The corresponding two-point 
PDFs of velocity differences were built on 2-4 x 10 9 pairs per 
snapshot each, depending on the pair separation. 

2.2. AMR Simulations 

The adaptive mesh refinement (AMR) simulation with ef- 
fective resolution of 2048 3 points was initialized by evolv- 
ing the flow on the root grid of 512 3 points for six dynamical 
times. Then one level of refinement by a factor of 4 was added 
that covers on average 50% of the domain volume. The grid 
is refined to better resolve strong shocks and to capture HD 
instabilities in the layers of strong shear. We use the native 
shock-detection algorithm of PPM, and we flag for refine- 
ment those zones that are associated with shocks with den- 
sity jumps in excess of 2. We also track shea r layers using 
the Fr obenius norm \\diUj(l - <5,y)||f, see also iKritsuk et al.1 
(2006). This second refinement criterion adds about 20% 
more flagged zones that would be left unrefined if only the 
refinement on shocks were used. The simulation was contin- 
ued with AMR for only 1 .2 td, which allows enough time for 
relaxation of the flow at high resolution, but does not allow 
us to perform time-averaging over many statistically indepen- 
dent snapshots. Therefore, we use the data from AMR simu- 
lations mostly to compare the quality of instantaneous statisti- 
cal quantities from simulations with adaptive and nonadaptive 
meshes. 

2.3. Random Forcing 

The initial velocity field is used, after an appropriate renor- 
malization, as a steady random force (acceleration) to keep 
the total kinetic energy within the box on an approximately 
constant level during the simulations. The force we applied 
is isotropic in terms of the total specific kinetic energy per 

dimension, u\ x = u\ y = u\ v but its solenoidal (V • uo.s = 0) 
and dilatational (V x uo.d = 0) components are anisotropic, 
since one of the three directions is dominated by the large- 
scale compressional modes, while the other two are mostly 
solenoidal. The distribution of the total specific kinetic en- 
ergy 



is defined as 



h = u • V x u. 



(4) 



£ =- I u z dV = £ s + £ D 



(3) 



between the solenoidal £$ and dilatational £q components is 
such that xo = £o.s/£o ~ 0.6. The forcing field is helical, but 

the mean helicity is very low, Iiq 2 <C h^, where the helicity h 



Note that in compressible flows with an isothermal equation 
of state the mean helicity is conserved, as in the incompress- 
ible case, since th e Ertel's potential vorticity is identically 
zero (lGaffetlll985l) . 

3. RESULTS 

In this Section we start with a general quantitative descrip- 
tion of our simulated turbulent flows and then derive their 
time-average statistical properties. We end up with assem- 
bling the pieces of this statistical picture in a context of 
a simple cascade model that extends the Kolmogorov phe- 
nomenology of incompressible turbulence into the compress- 
ible regime. 

3.1. Time-evolution of Global Variables 

The time variations of the rms Mach number and of the 
maximum gas density in the 1024 3 simulation are shown in 
Figs. [T]and |2] The kinetic energy oscillates between 18 and 
22, roughly following the Mach number evolution. Note the 
highly intermittent bursts of activity in the plot of p max {t). The 
time-average enstrophy 



and the Taylor scale 



A = 



1 



\u)\ 2 dV « 10 5 



;0.03 = 30A, 



(5) 



(6) 



where A is the linear grid spacing and u> = V x u is the vor- 
ticity. The rms helicity grows by a factor of 7.7 in the initial 
phase of the simulation and then remains roughly constant at a 
level of 1 .2 x 10 3 . While the conservation of the mean helicity 
is not built directly into the numerical method, it is still sat- 
isfied reasonably well. The value of h(t) is contained within 
±2% of its rms value during the whole simulation. 

If, instead of the Euler equations, we were to consider the 
Navier-Stokes equations with the explicit viscous terms, we 
could estimate the total viscous dissipation rate within the 
computational domain, 



where 



and 



e= $iUidV, 
h 



^Te^-¥-^ 



_ 1 / duj ^ duj 
IJ 2 \dxj dxi 



(7) 
(8) 
(9) 



is the symmetric rate-of-strain tensor, and the integration is 
done over the volume of the domain V = 1. Using vector 
identities, $ can be rewritten through the vorticity u) and the 
dilatation V • u as 



1 / 4 
* = Vxw--VV-B 

Re V 3 

and, by partial integration in equation (|7), 



€ = -Te (<|Vx«| 2 ) + l(|V.«| 2 ) 



(10) 



(11) 
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FIG. 3. — Probability density function for the gas density and the best-fit 
lognormal approximation. Note the excellent fit quality over eight decades 
in the probability along the high-density wing of the PDF. The sample size is 
2x 10". 

where the two terms on the right-hand side of equation ( fTTT i 
describe the mean dissipation rate due to solenoidal and di- 
latational velocities, respectively, e = e$ + co- 
in our Euler simulations, the role of viscous dissipation is 
played by the numerical diffusivity of PPM. This numerical 
dissipation does not necessarily have the same or similar func- 
tional form as the one used in the Navier-Stokes equations. 
It is still instructive, however, to get a flavor of the numeri- 
cal dissipation rate by studying the properties of the vorticity 
and dilatation. The so-calle d small-scale compressive ratio 
dKida & Orszagl[l990[[T992h . 



fcs : 



<|V-«| 2 ) 



l 2 ) + (|Vx«P) 



(12) 



represents the relative importance of the dilatational com- 
ponent at small scales. The time-average 7c$ sa 0.28. The 
mean fraction of the dissipation rate that depends solely on 
the solenoidal velocity component, es/e sw 0.65. Even though 
the two-dimensional shock fronts are dominating the geom- 
etry of the d ensity distribution in the dissipation range (see 
Section 13.81 ). the dissipation rate itself is dominated by the 
solenoidal velocity component that tracks corrugated shocks 
due to vortex stretching in the associated strong shear flows. 

3.2. Density PDF 

In isothermal turbulence the gas density does not corre- 
late with the local Mach number. As a result, the density 
PDF follows a lognormal distribution (rVaz auez-Semadeni 



1994tlPadoan et alfl 997: Passot & Vazquez-Semadeni ll998f 
Nordlund & Padoanll 19991 : lBiskampll2003l) . Fig. [3] shows our 
results for the time-average density PDF and its best-fit log- 
normal representation 



p(lnp)dlnp = 



\/2tt<j : 



x exp 



1 flnp-lnp 

~2 



dlnp, 



(13) 

where the mean of the logarithm of the density, In p, is deter- 
mined by 



lnp = -a z /2. 



(14) 



The lack of a Mach number-density correlation is illustrated 
in Fig. [4] where we show the two-dimensional PDF for the 
density and Mach number. The density distribution is very 
broad due to the very high degree of compressibility under 




log p 

FIG. 4. — Two-dimensional probability density function for the gas density 
and Mach number. The data here are based on a subvolume of 700 X 500 X 
250 ~ 10 s point s fro m a re presentative snapshot at t = 7 tj that we will also 
use later in Figs. 1 14l and ll7l The contours show the probability density levels 
separated by factors of 2. 

isothermal supersonic conditions. The density contrast ~ 10 6 
is orders of mag nitude higher than in the transonic case at 7 = 
1 .4 studied by lPorter etal] (120021) . Note the excellent quality 
of the lognormal fit in Fig. [3] particularly at high densities, 
over more than eight decades in the probability and a very 
low noise level in the data. 

If we express the standard deviation a as a function of Mach 
number M. as 

cr 2 =ln(l+b 2 M 2 ), (15) 

we get the best-fit value of b » 0.260 ±0.001 for log 10 p g 
J-2,2], which is sm aller than the b « 0.5 determined in 
Pad oan et all d 1997b for supersonic MHD turbulence. The 
powerful intermittent bursts in p max {t) (see Fig. [3]) correspond 
to large departures from the time-average PDF caused by 
head-on collisions of strong shocks. These events are usually 
followed by strong rarefactions that are seen as large oscilla- 
tions in the low-density wing of the PDF and also in the den- 
sity power spectrum. Intermittency is apparently very strong 
in supersonic turbulence. 

3.3. Velocity Power Spectra, Bottleneck Phenomenon, 
Numerical Dissipation, and Convergence 

We define the velocity power spectrum £(k) (or the specific 
kinetic energy spectral density) in terms of the Fourier trans- 
form of the velocity u 



Uik): 



1 



(2lT? 



u(x)e- 27Tik - x dx 



(16) 



as the square of the Fourier coefficient 

£(k)=^\u(k)\ 2 , (17) 

and then we define the three-dimensional velocity power 
spectrum as 

£(k)= [_£(k)8(\k\-k)dk, (18) 

where S(k) is the Dirac <5-function. The integration of the ve- 
locity power spectrum gives the specific kinetic energy (cf. 
eq. [3]), 

(19) 



£ = J £(k)dk=^(u 2 ) . 
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FIG. 5 . — Time-averaged velocity power spectra as a function of PPM dif- 
fusion coefficient K = 0.0, 0.1, and 0.2 at resolution of 512 3 grid points. The 
straight lines show the best linear fits to the spectra obtained for \ogk/k m j„ £ 
[0.6, 1.3]. The inertial range is barely resolved, since the driving scale over- 
laps with the bottleneck-contaminated interval. The slope of the "flat" part of 
the spectrum is mainly controlled by numerical diffusion. 

Our main focus in this Section is on the self-similar scaling 
of the power spectrum in the inertial range 

£(k) ~ k~ (20) 

that is limited by the kinetic energy input from the random 
force on large scales (kjk mm < 2) and by the spectrum flatten- 
ing in the near-dissipation part of the inertial range due to the 
so-called bottleneck effect related to a three-dimensional non- 
local mechanism of energy transfe r between modes of differ- 
ing length scales dFalkovichll 1 9941) . 

The bottleneck phenomenon has been observed both 
experimentally and in numerical simulations (e.g. , 
Porter etal.1 119941; iKaneda et all 120031: iDobler et all 120031: 
Haugen & B randenburg 2004). The strength of the bottle- 
neck depends on the way the dissipation scales with the 
wavenumber, and the effect is mo re pronounced w hen the 
dissipation grows faster than ~ k 2 ( Fa lkovichll 1 9941) . When 
we numerically integrate the Euler equations using PPM, the 
dissipation is solely det ermined by the method a nd affects 
scales smaller than 32 A ([Porter & Wood wardI 1 9941) . Accord- 
ing to iPorter etaf] (Il992al) . the effective numerical viscosity 
of PPM has a wavenumber dependency intermediate between 
~ k 4 and ~ k 5 . The high-order basic reconstruction scheme 
of PPM is designed to improve resolution in shocks and 
contact discontinuities and, therefore, numerical diffusion is 
controlled by various switches and is nonuniform in space. 

An addition of small diffusive flux near shocks 5 changes 
the scaling properties of numerical dissipation and reduces the 
flattening of the spectrum due to the bottleneck phenomenon. 
When the grid is not large enough to resolve the basic flow, the 
bottleneck bump in the spectrum is smeared and can be eas- 
ily misinterpreted as a shallower spectrum. Figure [5] shows 
how the slope of the "flat" section in the velocity power spec- 
trum at 512 3 depends on the value of the diffusion coeffi- 
cient K. The inertial range is unresolved in all three cases 
shown. The measured slope (3(K) changes from 1.91 ±0.01 
to 1.80 ±0.01 as the diffusion coefficient decreases from 0.2 
to zero and the bottleneck bump gets stronger. While the sta- 
tistical uncertainty of estimated power indices based on the 

s This is controlled by the parameter K in equation (4.5) in 
IColella & Woodward (19841) . 
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FIG. 6. — Compensated velocity power spectra from uniform grid PPM 
simulations at resolutions 256 3 , 512 3 , 1024 3 , and from AMR simulation with 
effective resolution of 2048 3 grid points. The wavenumbers are normalized 
to match k/kmin of the 1024 3 simulations at the Nyquist frequency. The dif- 
fusion coefficient K = for the uniform grid simulations, and K = 0.1 for 
the AMR simulation. All power spectra are time-averaged over t £ [6, 10] tj 
with the exception of the AMR one that is taken at t = 7.2 tj. The spectra 
demonstrate convergence for the inertial range of scales. 

fitting procedure is rather small, the variation of (3 as a func- 
tion of the diffusion coefficient is substantial and accounts for 
~ 33% of the difference between Burgers and Kolmogorov 
slopes! Clearly, it is difficult to get reliable estimates for the 
inertial range power indices from simulations with resolution 
up to 512 3 , since the slope of the spectrum is so dependent 
on numerical diffusivity. To study higher Mach number flows 
with PPM, an even higher resolution would be required. Note 
that with other numerical methods, one typically needs to use 
even larger grids, since the amount of dissipation provided 
by PPM is quite small compared to what is usually given in 
finite-difference schemes. In low-resolution simulations with 
high numerical dissipation that depends on the wavenumber 
as ~ k 2 or weaker, the power spectra appear steeper than they 
would have been when properly resolved. 

In simulations with AMR, the grid resolution is nonuni- 
form, and scaledependence (as well as nonuniformity) of nu- 
merical diffusion becomes even more complex. One should 
generally expect a more extended range of wavenumbers af- 
fected by numerical dissipation in AMR simulations than in 
uniform grid simulations, when the effective resolution is the 
same. The dependence of the effective numerical diffusiv- 
ity on the wavenumber in the AMR runs is, therefore, weaker 
than in simulations on uniform grids, and the bottleneck bump 
should be suppressed. 

In Figure [6] we combined the velocity power spectra from 
the unigrid simulations at 256 3 , 512 3 , and 1024 3 grid points 
with the diffusion coefficient K = to illustrate the conver- 
gence of the inertial range scaling with the improved reso- 
lution. All three simulations were performed with the same 
large-scale driving force, and the spectra were averaged over 
the same time interval, so the only difference between them is 
the resolution-controlled effective Reynolds number. We plot- 
ted the compensated spectra k 2 £(k) to specifically exaggerate 
the relatively small changes in the slopes. The over-plotted 
power spectrum from a single snapshot from our largest to 
date AMR simulation with the effective resolution of 2048 3 
grid points seems to confirm the (self-) convergence. Thus, 
one may hope that our estimates of the scaling exponents 
based on the time-averaged statistics from the 1024 simu- 
lation will not be too far off from those that correspond to 
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log k/k min 

FIG. 7. — Velocity power spectrum compensated by k 2 . Note a large-scale 
excess of power at £ £ [256, 1024]A due to external forcing, a short straight 
section in the uncontaminated inertial subrange I £ [40, 25 6] A, and a small- 
scale excess at i < 40A due to the bottleneck phenomenon. The straight 
lines represent the least-squares fits to the data for log k/k„,j„ £ [0.6, 1.1] and 
logVfcninS [1.2,1.7]. 
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FIG. 8. — Same as Fig. [7] but for the solenoidal velocity component. 
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FIG. 9. — Same as Fig. [7] but for the dilatational velocity component. 



Re — > oo. We further validate this statement in Section [3~4l 
where we discuss the scaling properties of the velocity struc- 
ture functions. Note also the suppression of the bottleneck 
bump visible in the AMR spectrum that occurs due to the 
higher effective diffusivity as we discussed earlier. 
Our reference velocity power spectrum from the 1024 3 sim- 



ulation is repeated in Fig.|7]with the best-fit linear approxima- 
tions. It follows a power law with an index (3 = 1.95 ±0.02 
within the range of scales £ € [40, 256] A and has a shal- 
lower slope, Pb = 1.72, for the bottleneck bump. After 
Helmholtz decomposition, £{k) = £s(k)+£o(k), the spectra for 
the solenoidal and dilatational components show the inertial 
range power indices of (3$ = 1.92 ±0.02 and (3d = 2.02 ±0.02, 
respectively, see Figs.[8]and|9] The difference in the slopes of 
£s(k) and £o(k) is about 5 a, i.e. significant. Both spectra dis- 
play flattening due to the bottleneck with indices of @b,s = 1 -70 
and /3b jj = 1 -74 in the near-dissipative range. The fraction of 
energy in dilatational modes quickly drops from about 50% at 
k = k m j„ down to 30% at kjk mm w 50 and then returns back to 
a level of 45% at the Nyquist frequency. 

In summary, the inertial range scaling of the velocity power 
spectrum in highly compressible turbulence tends to be closer 
to the Burgers scaling with a power index of (3 = 2 rather than 
to the Kolmogorov /? = 5/3 scaling sug gested for the mildly 
compressible transonic simulations by iPorter et al.1 |2002). 
The inertial range scaling exponents of the power spectra for 
the solenoidal and dilatational velocities are not the same, 
with the latter demonstrating a steeper Burgers scaling, con- 
taining about 30% of the total specific kinetic energy, and be- 
ing responsible for up to 35% of the global dissipation rate 
(see Section [3~T1 >. 

3.4. Velocity Structure Functions 

To substantiate this result, we studied th e scaling properties 
of the velocity structure functions (e.g.. iMonin & Yaglom 
[1971 

S p (£)=(\u(r+£)-u(r)\ p ) (21) 

of orders p = 1, 2, and 3, where the averaging (. . .) is taken 
over all positions r and all orientations of £ within the com- 
putational domain. Both longitudinal (u \\ £) and transverse 
(« _L £) structure functions can be well approximated by power 
laws in the inertial range 

S P (l)^f'\ (22) 
see Figures [T0lfl2l The low-order structure functions are less 
susceptible to the bottleneck contamination and m ight be bet- 
ter su ited for deriving the scaling exponents (Dobl eret al.1 
I20031) . 6 Let us now check how close the exponents derived 
from the data are to the K41 prediction ( p = p/3. 

It is natural to begin with the second-order structure func- 
tions, since the velocity power spectrum £(k) is the Fourier 
transform of S2CQ an d therefore (3 = £2 + 1 • The best-fit 
second-order exponents measured for the range of scales I be- 
tween 32A and 256 A, (J = 0.952 ± 0.004 and & = 0.977 ± 
0.008, are substantially larger than the K41 value of 2/3 and, 
within the statistical uncertainty, agree with our measured ve- 
locity power index (3 = 1 .95 ± 0.02. In homogeneous isotropic 

incompressible turbulence, is uniquely determined by s]J 
(and vice versa) via the first de Ka rman & Howarthl(ll938l) re- 
lation, 

Si=s\+ 1 -St (23) 
For the K41 inertial range velocity scaling, this translates into 
Si = \s\. (24) 

6 N ote, however, th at the bottleneck corrections grow with the or- 
der p IFalkovicr] fl994) and influence the structure functio ns in a nonlo- 
cal fashion, m ixing small- and large-scale information i Dobler et al] |2003l ; 
Davidson & Pearson 20Q5|). 
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FIG. 10. — First-order velocity structure functions and the best-fit power 
laws forlog 10 (^/A) e [1.5,2.5]. 
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FIG. 1 1 . — Same as Fig. 1101 but for the second-order structure functions. 
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Same as Fig. 1 1 01 but for the third-order structure functions. 



At Mach 6, the 2nd-order longitudinal and transverse struc- 
ture functions are also approximately parallel to each other in 

the (logS, logf)-plane, but the offset between them, S^/sf « 
1.27, is somewhat smaller than the K41 value of 4/3. This 
suggests that in highly compressible regimes some gener- 
alization of the first Karman-Howarth relation still remains 
valid, even though a substantial fraction of the specific kinetic 
energy (up to 1 /3 at Ai = 6) is concentrated in the dilatational 
component (cf. lMovallll95UISamtanev et al.ll200l . 




1 1.5 2 

FIG. 13. — Time-average total kinetic energy spectrum E(k) and solenoidal 
kinetic energy spectrum Es(k), both compensated by k 3 / 2 . The straight lines 
represent the least-squares fits to the data for \ogk/k„,j„ £ [0.5,1.2] and 
logk/k m i„ S [1.2, 1.8]. The inertial range slopes of the total and the solenoidal 
kinetic energy spectrum are the same within the statistical errors, which is 
not the case for the velocity power spectra. The dilatational component (not 
shown) also has the same slope. 



The first-order exponents Q = 0.533 ± 0.002 and = 
0.550 ±0.004 are again significantly larger than the K41 in- 
dex of 1/3 but also considera bly smaller than the index of 
1 for the Burgers model (e.g.. iFrisch & Becll2001h . Obser- 
vationally determined exponents for the molecular cloud tur- 
bulence are norm ally found within the range from 0.4 to 0.8 
(Brunt et al. 2003), in reasonable agreement with our result. 

We find that the third-order scaling exponents (J = 1 .26 ± 
0.01 and = 1.29 ±0.01 are also noticeably off from (3 = 1 
predicted by K41 in the incompressible limit. Our measure- 
ments for the three low-order exponents roughly agree with 
the estimates C j ~ ~ 0- 5, £2 ~ 0.9, and £3 w 1.3 obtained by 
iBoldvrev et al.l d2002l) from numerical simulations of isother- 
mal Mach 10 MHD turbulence at a resolution of 500 3 points. 
A rigorous result Cs = 1 holds for inc ompressible Navier- 
Stokes turbulence jKolmogorovlll941d, "4/5" law ) and for 
incompressible MHD dPolitano. & PouquetlfT998h . Strictly 
speaking, (3 = 1 is proved only for the longitudinal struc- 
ture function, S3 , (Kolm olForovll94icl) and for certain mixed 
structure functions (Politano, & Pouquetl ll998h . and in both 
cases the absolute value of the velocity difference is not taken 
(cf. eq. I2TI ). but still we know from numerical simulations 
that £3 « 1 in a compressible transonic re gime for the absolu te 
velocity differences, as in equation (f2TT> (IPorter etaDl2002). 

The fact that in supersonic turbulence £3 1 should be ac- 
counted for in extensions of incompressible turbulence mod- 
els into the strongly compressible regime. One approach is 
to consider relat ive (to the third order) scaling exponents as 
more universal (Dubrulle1 ll994h and refor mulate the mod- 
els in terms of the relative exponents (cf. [Boldvrevl 120021 : 
Bold vrev et alll2002l) . In many low-resolution simulations, 
exploiting the so-called extended self-similarity hypothesis 
dBenzi. et alj|1993l [l996h . the fo undations of which are not 
yet fu lly understood (see, e.g., ISreenivasan & B ershadskii 
2005), is the only resort to at least measure the relative ex- 
ponents. Another approach is to consider statistics of mixed 
quantities, like p v u, in place of pure velocity statistics naively 
inherited from incompressible regime studies. We shall dis- 
cuss the options for the former approach below and then re- 
turn to the latter in Sections [3.51 and [3791 
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The relative scaling exponents Z p = Cp/O of the first and 
the second order we obtain from our simulations are still 
somewhat higher than their incompressible nonintermittent 



analogues 1/3 and 2/3; z| - ^, 
0.76. This situation is similar to a small discrepancy be- 
tween the K41 predictions and expe rimental mea s ureme nts 
that stimulated Obu khovl (1 1962b and iKolmogorovl d 19621) to 
supplement the K41 theory by "intermittency corrections." If 
we follow a similar path in our strongly compr essible case and 
estima te the "corrected" exponents u sing the iBoldvrey et al.l 
(2002?) formula, which is similar to the lShe & Levequd(ll994l) 
model, but assumes that the most singular dissipative struc- 
tures are not one-dimensional vortex filaments but rather two- 
dimensional shock fronts, 



: 0.42, Zt = 0.43, and Z, 



- + 1- 



p/3 



(25) 



we get the exponents Z\ = 0.42 and Z2 = 0.74 that are very 
close to our estimates. One can argue though that at high 
Mach numbers, the dilatational velocities constitute a substan- 
tial part of the total specific kinetic energy, and therefore, a 
double set of singular velocity structures and their associated 
dynamics would lead to a compound Poisson statistic instead 
of a single Poisson process incorporated in equation $25[ (see 
She & Waymire 1995). To check whether this is indeed true 
or not, one needs to collect information about higher order 
statistics from simulations at resolution better than 1024 3 . 

3.5. Kinetic Energy Spectrum 

To explore the second option outlined in Section 3.4, we 
considered the energ y-spectrum function E(k) built on w = 
y/pu (e.g.. lLelelll994l) in addition to £(k) that depends solely 
on the velocity u and is traditionally used for characteriz- 
ing incompressible turbulence. We defined the kinetic energy 
spectral density E(k) as the square of the Fourier coefficient 



E(k) 



\\m\ 2 



determined by the Fourier transform 
1 



w(k) = 



w(x)e- 27rikx dx. 



(26) 



(27) 



(2tt) 3 Jv 

Then we used the standard definition of the three-dimensional 
kinetic energy spectrum function 



E(k)= / E(k)6(\k\-k)dk, 



(28) 



where S(k) as before is the <5-function. To avoid a possible 
terminological confusion, in the following we always refer to 
£(k) as the velocity power spectrum and to E(k) as the kinetic 
energy spectrum 



E(k)dk -- 



(29) 



At low Mach numbers, the scaling exponents of E(k) and 
£ (k) are nearly the same. However, while £(k) does not 
change much as the Mach number approaches unity (e.g., 
iPorter et"al . 2002), t he slope of E (k) gets shallower already 
by M. = 0.9 dKida & Orszadfl992l based on simulations with 
64 3 collocation points). Since at low Mach numbers the 
spectra are dominated by the solenoidal components, and the 
solenoidal velocity power spectrum hardly depends on the 




FIG. 14. — Enstrophy (|V X u| 2 ; top), density (middle), and "denstrophy" 
(|V X (y/~pu)\ 2 /p; bottom) distributions in a slice through the center of the 
subvolume of the computational domain at t = 7tj . The logarithmic gray- 
scale ramp is used to show the highest values in black and the lowest values 
in white. 

Mach number, the departure of the slope of E(k) from the 
slope of £(k) is mostly due to the sensitivity of the solenoidal 
part of the kinetic energy spectrum to the density-vorticity 
correlation that is felt already around M. ~ 0.5. Shallow 
slopes for E(k) with power indices around -3/2 were also 
detected in MHD si mulations of highly supersonic super- 
Alfvenic turbulence dLi et al.ll20"04l 512 3 grid points). 

In the inertial range, the time-averaged kinetic energy spec- 
trum scales as E(k) ~ k~ & with 0= 1.52 ±0.01 (see Fig.Qj). 
The spectrum displays flattening at high wavenumbers due to 
the bottleneck effect, which is similar to what we have already 



SIMULATIONS OF SUPERSONIC TURBULENCE 



9 




0.5 1 1.5 2 2.5 3 
'0910 ^Knin 

FIG. 15. — Time-average density power spectrum compensated by k. 
The straight lines represent the least-squares fits to the data for \ogk/k mm £ 
[0.6, 1.3] and \ogk/k min 6 [1.4, 1.7]. 

seen in the velocity power spectra £{k). The flat part scales 
roughly as ~ k~ 13 and occupies the same range of wavenum- 
bers as in the velocity spectrum. Since the kinetic energy 
spectrum is sensitive to the rather sporadic activity of den- 
sity fluctuations, time-averaging is essential to get a robust 
estimate for the power index. 

We performed a decomposition of w = ^fpu into the 
solenoidal and dilatational parts ws,d such that V-Wj = Vx 
Wd = 0, and computed the energy spectra for both compo- 
nents, E(k) = Es(k) + Eo(k). The inertial range spectral expo- 
nents for the solenoidal and dilatational components are the 
same, (3 = 1.52 ±0.01 for both, and also coincide with the 
slope of E(k). This remarkable property of E(k) suggests that 
in the supersonic regime 7 we are dealing with a single com- 
pressible cascade of kinetic energy where the density fluctua- 
tions provide a tight coupling between the solenoidal and di- 
latational mo des of w. This picture is s imilar in spirit to that 
discussed bv lKornrei ch & Scalo (2000, in t heir S ection 3.7). 
However, in contrast to iKornreich & Sca lo (2000), our simu- 
lations demonstrate that nonlinear HD instabilities are heavily 
involved in the kinetic energy transfer through the hierarchy 
of scales. The association of w with the kinetic energy dis- 
tribution as a function of scale makes this quantity a better 
candidate than the pure velocity u for employing in compress- 
ible cascade models, since w uniquely represents the physics 
of the cascade and there is no need to track the variation of 
Ed(K)/£(K) as a function of wavenumber in the inertial sub- 
range in addition to variations with the rms Mach number. 

The total kinetic energy power is dominated by the 
solenoidal component w$ over the whole spectrum. Within 
the inertial range, Es(k) contributes about 68% of the total 
power, then ~ 66% in the bottleneck-contaminated interval, 
and up to ~ 74% further down at the Nyquist frequency. Com- 
pare with the solenoidal part of the velocity power, £s(k), that 
constitutes about 65%-70% within the inertial range and only 
55% of the total power at the Nyquist frequency (see Section 

7 In our Mach 6 simulations at 1024 3 , the sonic scale £ s such that u(£ s ) ~ 
c = 1 is located in the middle of the bottleneck bump, and thus, the velocities 
within the inertial range are supersonic. One could possibly detect a break in 
the velocity power spectrum from a steep Burgers-like supersonic scaling at 
low k to a shallow Kolmogorov-like transonic scaling at higher wavenumbers 
in Mach 6 simulations with resolution of 4000 3 grid points or higher. Note, 
that as we show later in Section l3~9l the power spectrum of p'/ 3 u would not 
show such a break and would instead approximately follow the Kolmogorov 
5/3 law all over the inertial interval. 



1331) . 

Figure [14] illustrates the difference in structures seen in 
the enstrophy field {top) and in the "denstrophy" field (|V x 
(y/pu)\ 2 /p, bottom) that helps to understand the small-scale 
excess of power in E$(k) with respect to £s(k). While the cor- 
rugated shock surfaces (U-shapes), which are also the regions 
of very strong shear, are seen as dark wormlike structures of 
excess enstrophy and denstrophy in both the top and the bot- 
tom panels of Fig. [14] nearly planar shock fronts that carry 
a negligible amount of shear (this is also why they remain 
planar) are clearly missing in the enstrophy plot. Since the 
denstrophy field contains a greater number of sharper small- 
scale structures, it should be expected that E$(k) carries more 
small-scale power than £s{k). Overall, the structures captured 
as intense by the denstrophy field closely follow the regions 
of high energy dissipation rate (see the integrand in eq. JTJ 
and Fig. [17] bottom left). 

3.6. Density Power Spectrum 

The power spectrum of the gas density shows a short 
straight section with a slope of -1.07 ±0.01 in the range of 
scales from 250A down to 40A followed by flattening due 
to a power pileup at higher wavenumbers, see Fig. [15] Sim- 
ilar power-law sections, and excess of power in the same 
wavenumber ranges were seen earlier in the velocity and the 
kinetic energy power spectra, Figs. 171 and [T3l 

At the resolution of 512 3 , the bottleneck bump at high 
wavenumbers and the external forcing at low k leave essen- 
tially no room for the uncontaminated inertial range in k- 
space, even though the density spectrum at 512 3 (not shown) 
also has a straight power-law section. At 5 12 3 , the slope of the 
density power spectrum, -0.90 ± 0.02, is substantially shal- 
lower than -1.07 at 1024 3 . Thus, spectral index estimates 
based on low-resolution simulations bear large uncertainties 
due to the bottleneck contamination. We also note that the 
time-averaging over many snapshots is essential to get the cor- 
rect slope for the density power spectrum, since the density 
exhibits strong variations on very short (compared to tj) time- 
scales. The spectrum tends to get shallower after collisions of 
strong shocks, when the PDF's high-density wing rises above 
its average lognormal representation. 

There are good reasons to believe that in weakly com- 
pressible isothermal flows the three-dimensional density 
power spectrum s cales in the inertial subrange as ~ £~ 7 / 3 
ilBayly etaflll992l) . Our 512 3 transonic simulation at M = 1 
with purely solenoidal driving (not shown) returned a time- 
average power index of -1.6 ±0.1, which is probably still too 
shallow due to the unresolved inertial range (see also Sec- 
tion [3J]l^AnJndex of —1.7 at A4 = 1.2 was obtained by 
iKim & Rvul (|2005) at the same resolution (but with a more 
diffusive solver). Our simulations at higher resolution for 
M. = 6 give a slope of approximately —1.1. At even higher 
Mach numbers, the spectra tend to flatten further. From the 
continuity arguments, there should exist a Mach number value 
A^4i such that the power index of the density spectrum is ex- 
actly -5 /3, i.e. coincides with the Kolmogorov velocity spec- 
trum power index. Since apparently .M41 < 1, subsonic (but 
not too subsonic) compressible turbulent flows should have 
near-Kolm ogorov scaling, even for the density power spec- 
trum (cf. [Armstron g et allll995h . 

3.7. Coherent Structures in Supersonic Turbulence 

Let us now look more carefully at the morphology of those 
structures in supersonic turbulence that determine its scaling 
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FIG. 16. — Logarithm of the projected density from a snapshot of the AMR simulation with effective grid resolution of 2048 3 zones. The standard linear 
gray-scale ramp shows the highest density peaks in white and the most underdense voids in black. 



properties. The very first look at the projected density field 
shown in Figure[T6lreveals a plethora of filamentary structures 
very similar to what one usually sees as cirrus clouds in the 
sky. This is roughly what one can directly observe, albeit with 
a somewhat poorer resolution, in the intensity maps of nearby 
molecular clouds and in the WC-type turbulent stellar winds 
from the Wolf-Rayet stars. 

If one zooms in on a subvolume and takes a slab that is 
5 times thinner than the whole domain, one will start see- 
ing the elements of that cloudy structure, namely, numer- 
ous nested U- and V-shaped bow shocks or "Mach cones" 
dKritsuk et al.l |2006). Figure [T7] shows an example of co- 
herent structures that form in supersonic flows as a result 
of development and saturation of linear and nonlinear insta- 



bilities in shear layers formed by shock interactions. 8 The 
instabilities assist in breaking the dense post-shock gas lay- 
ers into supersonically moving fragments. Due to their ex- 
tremely high density contrast with the surrounding medium, 
the fragments behave as effectively solid obstacles with re- 
spect to the medium. These coherent structures are transient 
and form in superson ic collisions of counter-pr o pagating flow 
patches or "blobs" feritsuk & Normanl 120041: iHeitsch et all 
l200l iKritsuket all [2006). Then they dissipate and similar 
structures appear again and again at other locations. The 

8 The (incomplete) list of potential candidates includes the Kelvin- 
Helmholtz instability (Helmholtz 1868), the nonlinear supersonic vortex 
sheet instability (Miles 1957), and the nonlinear thin shell instability 
IVishniad 1994ft . 
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FIG. 17. — Coherent structures in Mach 6 turbulence at resolution of 1024 3 . Projections along the minor axis of a subvolume of 700 X 500 X 250 zones for the 
density (top left), the enstrophy (|V X u| 2 ; top right), the dissipation rate (the integrand |$*if| in eq.[7j bottom left), and the dilatation (V ■«; bottom right). The 
logarithmic gray-scale ramp shows the lower values as dark in all cases except for the density. The inertial subrange structures correspond to scales between 40 
and 250 zones and represent a fractal with D,„ pb 2.4. The most singular structures in the dissipation range (I < 30A) are shocks with fractal dimension D m = 2. 



nested V-shapes can be best seen in slabs with a finite thick- 
ness of the order of their characteristic size (i.e. ~ 200A 
in our 1024 3 run) that belongs to the resolved inertial sub- 
range of scales. They are less noticeable in full-box projec- 
tions where multiple structures of the same sort tend to over- 
lap along the line of sight as in Figure [16] Since the "cones" 
are essentially three-dimensional, they can also be hardly seen 
in thin slices, see Fig. [14] Patterns of the same morphol- 
ogy were identified in the peripheral parts of the Ml -67 neb- 
ula observed in H a e mission by the Hubble Space Telescope 
dGrosdidieretal.1120011) . 

3.8. Fractal Dimension of the Mass Distribution 

It is known from observations and experiments that the dis- 
tribution of dissipation in turbulent flows is intrinsically in- 
termittent and this intermittency bears a hierarchical nature. 
This fact is in apparent contradiction with the K41 theory that 
assumes that the rate of dissipation is uniform in space and 
constant in time. Mande lbrot] d 19741) introduced a new con- 
cept of the "intrinsic fractional dimension of the carrier," D, 
that characterizes the geometric properties of a subset of the 
whole volume of turbulent flow where the bulk of intermittent 
dissipation occurs. He also studied the relation between D and 
scaling properties of turbulence. In particular, he suggested 
that isosurfaces of scalars (such as concentration or tempera- 
ture) in turbulent flows with Burgers and Kolmogorov statis- 



tics are best describe d by fractal dimen sions of 3- 1/2 and 
3-1/3, respectively (Mandelbrot 1975). In the following we 
focus on similar questions with respect to supersonic turbu- 
lence. 

In turbulence research there are several approaches to assess 
the flow geometry quantitatively. Indirectly, the information 
on fractional dimension of dissipative structures in turbulent 
flows can be obtained as a by-product of application of phe- 
nomenological cascade models, e.g. the hierarchical s truc- 
ture (HS) mod el developed by I She & Levequd (1 19941) and 
Dubrulle (1994) for incompressible Navier- Stokes turbulenc e 
and later extended to supersonic flows by Bold yrevl (120021) . 
The HS model is based on the assumption of log-Poisson 
statistics for the energy dissipation rate and provides an esti- 
mate for fractal dimension, D, of the most singular structures 
based on the detailed knowledge of high-order velocity struc- 
ture fu nctions (e.g., iKritsuk & Nor man 2004; Padoa rTet al.l 
2004b). While the HS model appears to reproduce experi- 
mental data quite well in a rather diverse set of applications 
extending far beyond the l imits of the o riginal study of fully 
developed turbulence (e.g. jLiu et alj2004h . the uncertainty in 
the scaling properties of high-order velocity statistics remains 
large. To nail down the D value thus obtained from numerical 
experiments, one has to go beyond the current resolution lim- 
its to resolve the inertial subrange dynamics sufficiently well 
and to provide a statistical data sample of a reasonably large 
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size. That would reduce the effects of the statistical noise 
and anisotropics which tend to strongly co ntaminate the high- 
order velocity statistics (e.g.. iPorter et al.|[2002l) . We leave a 
more detailed discussion of the HS model to a follow-up paper 
and focus here instead on a direct measurement of the fractal 
dimension of the density field using conventional techniques. 
Note, however, that there is no trivial physical connection 
between the fractal dimensions measured by these two tech- 
niques, i.e. they do not necessarily refer to the same fractal 
object since one is based exclusively on the velocity informa- 
tion, while the other is based on the density information. 

For physical systems in three-dimensional space, like su- 
personic turbulent boxes, it is often easier to directly estimate 
the mass dimension that contains information about the den- 
sity scaling with size. To do so, we select a sequence of box 
sizes l\ > ii > ... > £„ and cover high-density peaks with con- 
centric cubes of those linear sizes. Denoting M{£/) as the mass 
contained inside the box of size one can plot logM(^,) ver- 
sus log 1; to reveal a scaling range where the points follow a 
straight line, 

M(£) oc £ D ™ . (30) 

The slope of the line, D m , then gives an estimate of the mass 
dimension for the density distribution. Note that the centers 
of the box hierarchies cannot be chosen arbitrarily. They must 
belong to the fractal set, or else, in the limit £ — ► 0, one would 
end up with boxes containing essentially no mass. In prac- 
tice, we selected as centers the positions where the gas den- 
sity is higher or equal to one-half of the density maximum for 
a given snapshot. We then determined D m by taking an aver- 
age over all centers from five statistically independent density 
snapshots separated from each other by roughly one dynami- 
cal time. 

Figure [18] shows the mass scaling with the box size I com- 
pensated by £ ~ 2 based on the procedure described above. The 
slope of the mass-scale dependence breaks from roughly 2 
in the dissipation-dominated range of small scales to ~ 2.4 
in the inertial subrange that is free from the bottleneck effect. 
The transition point at £ ~ 30A coincides with the break point 
present in all power spectra discussed above. On small scales, 
the geometry of the density field in supersonic turbulence is 
dominated by two-dimensional shock surfaces. This is where 
the numerical dissipation dominates in PPM. 

On larger scales, where the dissipation is practically negli- 
gible, the clustered shocks form a sophisticated pattern with 
fractal dimension D m = 2.39 ±0.01. Since the dynamic range 
of our simulations is limited, we can only reproduce this scal- 
ing in a quite narrow interval from ~ 40 A to ~ 160 A, and 
thus, the actual uncertainty of our estimate for D m in the in- 
ertial range can still be on the order of ±0.1 or even larger. 
The break present in the (logM-log^) relation may indicate 
a change in intermittency properties of turbulence at th e meet- 
ing po int of the inertial and dissipation ranges (see, e.g. lFrischl 
fl995h . 

3.9. A Simple Compressible Cascade Model 

More than half a century ago. Ivon Weizsackerl (1195 lb intro- 
duced a phenomenological model for three-dimensional com- 
pressible turbulence with an intermittent, scale-invariant hier- 
archy of density fluctuations described by a simple equation 
that relates the mass density at two successive levels to the 
corresponding scales through a universal measure of the de- 




0.5 1 1.5 2 2.5 

[og 10 e/zl 

FIG. 18. — Gas mass within a box of size I normalized by I 2 as a function 
of the box size. The mass dimension D„, of the density distribution is deter- 
mined by the slope of this relationship. A horizontal line would correspond 
to D m = 2. The data points represent averages over five density snapshots. 
The straight lines are the least-squares fits to the data for log£/A £ [0.5, 1.2] 
and log ijA e [1.7,2]. 

gree of compression, a, 

Pu-l \£u-lj 

The only free parameter of the model is the geometrical 
factor a which takes the value of 1 in a special case of 
isotropic compression in three dimensions, 1/3 for a perfect 
one-dimensional compression, and zero in the incompressible 
limit. 

The kinetic energy supplied to the system at large scales 
is being transferred throu gh the hierarchy by nonlinear inter- 
actions. LighthijJ (1 19551) pointed out that, in a compressible 
fluid, the mean volume energy transfer rate pu 2 u /£is constant 
in a statistical steady state, so that 

u~{i/p) x '\ (32) 

From t hese two equations, assuming mass conservation. lFleckl 
( 1996) derived a set of scaling relations for the velocity, spe- 
cific kinetic energy, density, and mass: 9 

m~^ 1/3+q , (33) 

£(k)~k- ~/<r 5/3 - 2Q , (34) 

P~r 3a , (35) 

M(£)~£ D » ~^ 3 " 3q , (36) 

where all the exponents depend on the compression measure 
a. The compression measure, in turn, is a function of the rms 
Mach number of the turbulent flow. In the incompressible 
limit, /A — ► and a — ► 0. There are also scaling relations for 
v = p l l^u that follow directly from equation (f32t and extend 
the incompressible K41 velocity scaling into the compressible 
regime, 10 

v p = (p xl \) p ~l pl \ (37) 

These hint at a unique generalization of the velocity structure 
functions for compressible flows, 

S p (£) = (\v(r+D-v(r)\ p ) ~ £"'\ (38) 

9 Similar ideas were discussed earlier by Biglari & Diamond 1 1988, 1989) 
in a slightly different context of self-gravitating compressible fluids. 

10 Note that Kolmogorovl <1941al lbT) and|Obukhov 1 194 J) only considered 
p < 3 and never went as far as to consider p > 3 (cf. Frisch 1995). 
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with S3 ~ £ (see also the discussion in Section l3~4t . The scal- 
ing laws expressed by equation (l38l should not necessarily be 
exact and, as the incompressible K41 scaling, may require 
"intermittency corrections" that will be addressed in detail 
elsewhere. Our discussion here is limited to low-order statis- 
tics (p < 3) for which the "corrections" are supposedly small. 
The compression measure a could also depend on the order p 
in addition to its Mach number dependency. However, as we 
show below, a linear approximation a = const built into this 
simple cascade model may prove reasonable for the low-order 
density and velocity statistics. 

Before turning to properties of the density-weighted veloci- 
ties V, let us first assess the predictive power of Fleck's model 
using the statistics we already derived. Since from the simula- 
tions we know the scaling of the first-order velocity structure 
functions, we can get an estimate of the geometric factor a for 
the Mach 6 flow using equation d33l . Assuming the inertial 
range scaling, Si ~ £ ()M , we get a w 0.21 and D m ss 2.38 from 
equation d36l ). This is consistent with our direct measurement 
of the mass dimension, D m sa 2.4, for the same range of scales. 

We can also do a consistency check based on the second- 
order statistics. From the inertial range scaling of the second- 
order structure functions, we get £(k) ~ k~ 197 and therefore 
a = 0.15. This value corresponds to the fractal dimension 
D m « 2.55 that is higher, but still reasonably close to our di- 
rect estimate. The discrepancy can in part be attributed to a 
deficiency of Fleck's model that asymptotically approaches 
the space-filling K41 cascade in the limit of weak compress- 
ibility (a — > 0) and thus does not take the full account of the 
intermittency of the velocity field. 

It is inter esting to compare the inte r face d imension D,- in- 
troduced by Meneveau & Sreenivasan ( 1990) for intermittent 
incompressible turbulence based on the so-called Reynolds 
number similarity 

A = 2 + 6, (39) 

where £i as before is the exponent of the first-order veloc- 
ity structure function, with the mass dimension D m = 3 -3a. 
Since both quantities characterize the same intermittent struc- 
tures, they may well refer to the same fractal object in our 
compressible case. If one assumes = D m = D, keeping in 
mind that £j = 1/3 + a as follows from d33l . it is easy to com- 
pute the compression parameter a = 1/6 ~ 0.16(6) as well as 
the fractal dimension, D = 2.5, t hat appear to be in accord with 
the original proposal by Fleck (1 19961) and with the numbers 
listed above. 

Finally, we can use the scaling exponents of the velocity 
power spectrum £(k) ~ k~ and of the kinetic energy spec- 
trum E(k) ~ k~ to find a from the density scaling in equa- 
tion J35b . We get 



E(k) 
£{k) 



.0.45 



-3a 



(40) 



and thus, a = 0. 15, consistent with the previous estimate based 
on the second-order statistics. 

Overall, within the uncertainties, Fleck's model appears 
to successfully reproduce the low-order velocity and density 
statistics from our numerical simulations and supports our di- 
rect measurement of the mass dimension D,„ as well. 

Let us now check how well equation d37l i is satisfied in 
our simulations by computing the power spectrum £(£) of 
the density-weighted velocity v as a p = 2 diagnostic. We 
define £(£) in exactly the same way as £(k) was defined in 
Section 13.31 but substitute v in place of the velocity u. As 



V 




1 1.5 

'0910 ^Knin 

FIG. 19. — Time-averaged power spectrum of the density-weighted velocity 
V = p'/ 3 H compensated by £ 5 / 3 . The straight lines represent the least-squares 
fits to the data for logifc/fe™ 6 [0.5,1.1] and \ogk/k mi „ 6 [1.2,1.8]. The 
inertial subrange slope is in excellent agreement with the model prediction. 



can be seen from Figure [19] the inertial range power index, 
j3 = 1 .69 ± 0.02, is very close to the Kolmogorov value of 5/3. 

To further check our conjecture expressed by equation [38] 
we measured the first three scaling exponents for the modi- 
fied structure functions S p using a subsample of 35 snapshots 
evenly distributed in time through t € [6, 10]f</ using a sam- 
ple of 2 x 10 9 point pairs per PDF per snapshot. The result- 
ing third-order exponents are very close to unity as expected, 

& = 1.01 ± 0.02 and (J = 0.95 ± 0.02. This strongly sug- 
gests that a relationship similar to Kolmogorov's "4/5" law 
may hold for compressible flows as well. 

Both the first-order exponents Ct = 0.467 ± 0.004 and (J 1 = 
0.451 ±0.004 and the second-order exponents = 0.80 ± 

0.01 and (J = 0.76 ±0.01 slightly deviate from the K41 p/3 
scaling as would be the case in intermittent turbulence. The 
corresponding relative (to the 3rd order) exponents = 0.46 

and zj 1 = 0.47 and = 0.79 and zJj = 0.80 are also some- 
what higher than their counterparts for the velocity structure 
functions discussed in Section 13.41 indicating certain struc- 
tural differences between the u and v fields. 

Note, as the Mach number changes from M. < 1 to M. = 6, 
the slope of the density power spectrum gets shallower from 
-7/3 to —1.1, but the slope of the velocity power spectrum 
gets steeper from -5/3 to —1.9. At the same time, the power 
spectrum of the mixed variable v remains approximately in- 
variant, as predicted by the model. 

Our numerical experiments thus confirm one of the ba- 
sic assumptions adopted in the compressible cascade model, 
namely, that equation ([32l for the energy transfer rate holds 
true. The first assumption concerning the properties of the 
self-similar hierarchical density s tructure, equation (l3"TT i put 
forward by von Weizsiicker (1951), seems to be also satisfied 
in the simulations quite well, at least to the first order. 

4. DISCUSSION 

The major deficiency of the numerical experiments dis- 
cussed above is still their limited spatial resolution that 
bounds the integral scale Reynolds numbers to values much 
smaller than those estimated for the real molecular clouds. 
This hurdle apparently cannot be overcome in the near future, 
but still the progress achieved in the past 15 years in this di- 
rection is very impressive. 



14 



KRITSUK, NORMAN, PADOAN, & WAGNER 



The second important deficiency of our model is the lack 
of magnetic effects which are known to be essential for star 
formation applications. This subject still remains a topic for 
the future work awaiting the development of a high-quality 
MHD solver suitable to modeling of supersonic flows at mod- 
erately high Reynolds numbers with computational resources 
available today. 

Another set of potential issues relates to the external driv- 
ing force that is supposed to simulate the energy input by 
HD instabilities in real molecular clouds. The large-scale 
driving force we used in these simulations is not perfectly 
isotropic due to the uneven distribution of power between 
the solenoidal and dilatational modes (perhaps, a typical sit- 
uation for the interstellar conditions). We also use a static 
driving force that could potentially cause some anomalies on 
timescales of many dynamical times. However, while strong 
anisotropie s can significantly affect the scaling o f high-order 
moments dPorter et al.ll20Q2i: iMininni et al.ll2006l) . the depar- 
tures from Kolmogorov-like scaling we observe in the lower 
order statistics appear to be too strong to be explained solely 
as a result of the specific properties of the driving. The sensi- 
tivity of our result to turbulence forcing remains to be verified 
with future high-resolution simulations involving a variety of 
driving options. 

The options for observational validation of our numerical 
models are limited. Interstellar turbulence in general and su- 
personic turbulence in molecular clouds in particular coul d 
hardly be uniform and/or isotropic ( Kap lan & Pikelnerl 19701) . 
There are multiple driving mechanisms of different natures 
operating on different scale s in the ISM dNorman & Ferraral 
1996; Mac Low & Klessenl 120041) . so the source function of 
turbulence is expected to be broadband. Various observational 
techniques employed to extract information about the scaling 
properties of turbulence have their own limitations, includ- 
ing a finite instrumental resolution, insufficiently large data 
sets, inability to fully access the three-dimensional informa- 
tion without additional a priori assumptions, etc. Moreover, 
complexity of the effective equation of state of the ISM and 
many other physical processes so far ignored in the simplified 
numerical models should also make the comparison of ob- 
servations and simulations uncertain. Nevertheless, it makes 
sense to compare our results with the scaling properties of su- 
personic turbulence obtained from observations. 

Applying the velocity channel a nalysis (VCA) technique 
(Lazarian & Pogosvan 2000, 2006) to power spectra of inte- 
grated intensity maps and single-vel ocity channel maps of the 
Perseus region, iPadoan et al.1 (120061) found a velocity power 
spectrum index (3 = 1.8 that is reasonably close to our mea- 
surement (3 = 1.95. The struc ture function exponents mea- 
sured for the Ml -67 nebula by iGrosdidier et alJ (1200 ll) . Ci ~ 
0.5 and (2 ~ -9, match quite nicely with our results discussed 
in SectionES Ci = 0.54 and C2 = 0.97. 

Using maps of the 13 CO / = 1 - emission line of the 
molecular cloud com plexes in Perseus, Taurus, and Rosetta, 
IPadoan et al.1 (12004 al) computed the power spectra of the col- 
umn density estimated using the LTE 11 method (Dickman 
1978). The slopes of the measured spectra corrected for tem- 
perature and saturation effects on the 13 CO 7=1-0 line, 
-0.74 ±0.07, -0.74 ±0.08, and -0.76 ±0.08, respectively, 
are notably shallower than our estimate for the density spec- 
trum power index -1 .07 ± 0.01 (see Section [3~6l ). The appar- 
ent 4 a discrepancy is most probably due to an insufficiently 

' 1 Local thermodynamic equilibrium. 



high Mach number adopted in our simulations. Alternatively, 
it can be attributed to anisotropies in the molecular cloud tur- 
bulence, intermittent large-scale driving force acting on the 
clouds, or limitations of the LTE method. 

As far as the fractal dimension is concerned, the 
observational measurements for molecular clouds and 
star- forming regions tend to cluster around D = 2.3 ± 
0.3 dE lmegreen & Falgarone 1996; Elmegreen & Elmegreen 
12001 " Stochastically variable turbulent [WC] winds from 
"dustars" feeding the ISM also demonstrate similar fractal di- 
mension D = 2.2-2.3, and the same morphology of clumps 
forming and dissipating in real time is clearly seen in the outer 
parts of the wind where the picture is not smeared by projec- 
tion effects (Grosdidi er et al.ll200l1) . 

5. CONCLUSIONS 

Using large-scale numerical simulations of nonmagnetic 
highly compressible driven turbulence at an rms Mach num- 
ber of 6, we were able to resolve the inertial range scaling and 
have demonstrated that: 

1. The probability density function of the gas density is 
perfectly represented by a lognormal distribution over 
many decades in probability as predicted from simple 
theoretical considerations. 

2. Low-order velocity statistics deviate substantially from 
Kolmogorov laws for incompressible turbulence. Both 
velocity power spectra and velocity structure functions 
show steeper than Kolmogorov slopes, with the scaling 
exponents of the third-order velocity structure functions 
far in excess of unity. 

3. The density power spectrum is instead substantially 
shallower than in weakly compressible turbulent flows. 

4. The kinetic energy power spectrum (built oiiwe ^fpu) 
is shallower than Kolmogorov's 5/3-law. The power 
spectra for both the solenoidal and the dilatational parts 
of w obtained through Helmholtz decomposition (w = 
Wj+Wd) have the same slopes as the total kinetic en- 
ergy spectrum pointing to the unique energy cascade 
captured by E(k); their shares in the total energy bal- 
ance are 68% and 32%, respectively. 

5. As should have been expected, the mean volume energy 
transfer rate in compressible turbulent flows, pu 2 u / 1, is 
very close to constant in a (statistical) steady state. Our 
simulations show that the power spectrum of v = p l ^u 
scales approximately as k~ 5 / 3 and the third order struc- 
ture function of v scales linearly with £ in the inertial 
range. 

6. The directly measured mass dimension of the "fractal" 
density distribution is about 2.4 in the inertial range and 
2.0 in the dissipation range. The geometry of super- 
sonic turbulence is dominated by clustered corrugated 
shock fronts. 

These results strongly suggest that the Kolmogorov laws 
originally derived for incompressible turbulence would also 
hold for highly compressible flows as soon as they are ex- 
tended by replacing the pure velocity statistics with statistics 
of mixed quantities, such as the density-weighted fluid veloc- 
ity V = p 1//3 M. 
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Both the scaling properties and the geometry of supersonic 
turbulence we find in our numerical experiments support a 
phenomenological description of intermittent energy cascade 
in a compressible turbulent fluid and suggest a reformulation 
of inertial range intermittency models for compressible turbu- 
lence in terms of the proposed extension to the K41 theory. 

The scaling exponents of the velocity diagnostics, the mor- 
phology of turbulent structures, and the fractal dimension 
of the mass distribution determined in our numerical ex- 
periments demonstrate good agreement with the correspond- 
ing observed quantities in supersonically turbulent molecular 
clouds and line radiation-driven winds from carbon-sequence 



Wolf-Rayet stars. 
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